	SUBROUTINE SHANGPU(X,N,M,K0,SE,T,OMEGA,SIGMA2,FAI,B)
	INTEGER::N,M
	REAL(8),DIMENSION(N)::X
	REAL(8),DIMENSION(N-1,N-1)::B
	REAL(8),DIMENSION(0:N-1)::SIGMA2
	REAL(8),DIMENSION(N-1)::FAI
	REAL(8),DIMENSION(0:M)::SE,T,OMEGA
	REAL(8)::A1,A2,A3,A4
	REAL(8)::PI
	PI=4*ATAN(1.)
	B=0
	A3=0
	DO I=2,N
		A3=A3+X(I)*X(I-1)
	END DO
	A4=(X(1)*X(1)+X(N)*X(N))/2
	DO I=2,N-1
		A4=A4+X(I)**2
	END DO
	B(1,1)=A3/A4
	DO K=1,N-2
		A1=0
		A2=0
		DO I=K+2,N
			A3=X(I)
			A4=X(I-K-1)
			DO J=1,K
				A3=A3-B(K,J)*X(I-J)
				A4=A4-B(K,J)*X(I-K-1+J)
			END DO
			A1=A1+A3*A4
			A2=A2+A3**2+A4**2
		END DO
		B(K+1,K+1)=2*A1/A2
		DO J=1,K
			B(K+1,J)=B(K,J)-B(K+1,K+1)*B(K,K+1-J)
		END DO
	END DO
	SIGMA2(0)=0
	DO I=1,N
		SIGMA2(0)=SIGMA2(0)+X(I)*X(I)
	END DO
	SIGMA2(0)=SIGMA2(0)/N
	DO K=1,N-1
		SIGMA2(K)=SIGMA2(K-1)*(1-B(K,K)**2)
		FAI(K)=SIGMA2(K)*(N+K)/(N-K)
	END DO
	DO K=3,N-1
		IF(FAI(K-1)>FAI(K).AND.FAI(K)<FAI(K+1))THEN
			K0=K
			GOTO 10
		END IF
	END DO
	K0=N-1
10	DO L=0,M
		A3=0
		A4=0
		DO K=1,K0
			A3=A3+B(K0,K)*COS(PI*L*K/M)
			A4=A4+B(K0,K)*SIN(PI*L*K/M)
		END DO
		SE(L)=SIGMA2(K0)/((1-A3)**2+A4**2)
		OMEGA(L)=PI*L/M
	END DO
	DO L=1,M
		T(L)=2.*M/L
	END DO
	END
! 7.5.6例
! 	以北京1951年—1980年1月平均温度为例。序列长度为N=30,取最大波数M=13。
	PROGRAM MAIN
	INTEGER,PARAMETER::N=26  	  !序列长度
	INTEGER,PARAMETER::M=13		  !最大波数
	REAL(8),DIMENSION(N)::X
	REAL(8),DIMENSION(0:M)::SE,T,OMEGA
	REAL(8),DIMENSION(N-1)::FAI
	REAL(8),DIMENSION(0:N)::SIGMA2
	REAL(8),DIMENSION(N-1,N-1)::B
	REAL(8)::XBAR
	OPEN(10,FILE='BEIJING.DAT')
	READ(10,*)X
	CLOSE(10)
	N0=N
	XBAR=0
	DO I=1,N
		XBAR=XBAR+X(I)
	END DO
	XBAR=XBAR/N
	X=X-XBAR
	CALL SHANGPU(X,N,M,SE,T,OMEGA,SIGMA2,K0,FAI,B)
	OPEN(12,FILE='SHANGPU.DAT')
	WRITE(12,'(" B=")')
	DO I=1,K0
		WRITE(12,'(<K0>E10.3)')(B(I,J),J=1,K0)
	END DO
	WRITE(12,'(" K0=",I4)')K0
	WRITE(12,'(3X,"I",7X,"SE",10X,"T",9X,"OMEGA")')
	DO I=0,M
		WRITE(12,'(1X,I3,3E12.4)')I,SE(I),T(I),OMEGA(I)
	END DO
	WRITE(12,'(2X," I",3X,"SIGMA2(I)",5X,"FAI(I)")')
	DO I=1,N-1
		WRITE(12,'(I4,2E12.4)')I,SIGMA2(I),FAI(I)
	END DO
	CLOSE(12)
	END
